#######做eQTL富集分析
case=read.table("bias_AShM_BF_no_motif.txt",head=T,sep="\t")
con=read.table("nobias_AShM_total_reads10.txt",head=T,sep="\t")

case=case[case$BF_in_DC>10,]
case_id=paste(case$unitID,case$Ref,case$Alt,sep="_")
caseid=gsub("chr","",case_id)
caseid=gsub(":","_",caseid)

conid=gsub(":","_",as.character(con$unitID))
conid=gsub("chr","",conid)

file=read.table("/mnt/md1200/5/zhaocunyou/wangzhongju/eQTL_analysis_by_wangzhongju/eQTL_data/GTEx_Analysis_v6p_all-associations/Prostate_Analysis.v6p.all_snpgene_pairs.txt",head=T,sep="\t")
file=tidyr::separate(file,variant_id,into=c("chrome","pos","ref","alt","what"),sep="_")
file$unitID=paste(file$chrome,file$pos,file$ref,file$alt,sep="_")

case=file[file$unitID %in% caseid,]
con=file[file$unitID %in% conid,]

length(conid)
[1] 12284

length(caseid)
[1] 727

> length(unique(case$unitID))
[1] 524
> length(unique(con$unitID))
[1] 9424
con_sig=con[con$pval_nominal<0.05,]
con_nosig=con[con$pval_nominal>0.5,]
con.nosig.id=unique(con_nosig$unitID)
con.sig.id=unique(con_sig$unitID)
 test1=setdiff(con.nosig.id,con.sig.id)
 
 case_sig=case[case$pval_nominal<0.05,]
case_nosig=case[case$pval_nominal>0.5,]
case.nosig.id=unique(case_nosig$unitID)
case.sig.id=unique(case_sig$unitID)
 test2=setdiff(case.nosig.id,case.sig.id)

fisher.test(matrix(c(length(case.sig.id),length(test2),length(con.sig.id),length(test1)),nrow=2))$estimate
fisher.test(matrix(c(length(case.sig.id),length(test2),length(con.sig.id),length(test1)),nrow=2))$p.value

415	727	7485	12284